Paige Duffin Data Analysis Project
Working title- will give proper title later
Paige Duffin
2019-10-11
Organizing this into tabs for the time being
Part 1: Introduction
Summary/Abstract Write a summary of your project.
General Background Information
Outbreaks of infectious disease are an important, yet poorly understood, driving force in population biology and, while prevalent in both terrestrial and marine systems, the manner by which outbreaks influence biotic relationships, age demographics, community structure and function, and trophic interactions in the aquatic realm consistently lags behind that of terrestrial disease ecology. The growing body of evidence in this field, however, supports the observation that marine disease epidemics are increasing in frequency and severity (Orth et al., 2006; Waycott et al., 2009). Elucidating the relationship between aquatic species and their pathogenic diseases is pertinent in the field of marine ecology because infection outbreaks have the potential to drastically alter ecosystem functionality (Blakesley et al., 2002; Burge et al., 2013). Many of the organisms which have faced chronic and/or severe outbreaks of disease (such as sea urchins, scleractinian corals and seagrasses) are also considered keystone species and/or ecosystem engineers, meaning that they contribute highly valuable services or functions to their surrounding community and ecosystem. Thus, it comes as no surprise that disease-driven mass mortalities of these species often generate waves of ecological permutation, ranging from temporary local disruptions to permanent phase shifts (Burge et al., 2013), such as the transformation from a coral to macroalgal dominated Caribbean reef structure that accompanied the massive Diadema antillarum (black sea urchin) die-off of the early 1980s (Lessios et al., 1984; Lessios, 1988). This event, likely caused by a biological pathogen (Schultz et al., 2016), reduced Diadema populations to less than 1% of their original size (Lessios, 1988). Mass mortality events impacting critical foundation species, ecosystem engineers, or keystone species such as the black sea urchin have been coined ‘marine disease emergencies’ due to the detrimental cascade of events that often succeeds them (Miner et al., 2018).
Sea Star Wasting Disease (SSWD) is another epizootic crisis facing modern day coastal ecosystems. Outbreaks of asteroid wasting have been documented periodically since the late 1970s and SSWD describes a suite of symptoms observed across a broad range of sea star species, most of which play integral parts in shaping their community structure. Generally speaking, the wasting disease events which punctuated the past four decades were relatively brief, in localized areas, and largely failed to capture the attention of the scientific community. Beginning in summer 2013, however, mass mortalities of sea stars due to wasting disease have caused unprecedented damage, owing to the geographical and temporal extent of impact. This ongoing epizootic, which has killed millions of asteroids across over 20 taxa, is widely referred to as the largest disease event sweeping through a wildlife marine species in documented history (Hewson et al. 2014, Gudenkauf & Hewson, 2015, Eisenlord et al. 2016).
Description of data and data source
Describe what the data is, what it contains, where it is from, etc.
The Source
The vast majority of disease documentation and surveys cover the west coast of the United States. While it is almost certainly true that SSWD is most prevalent in this region, there is a strong bias towards the discovery of diseased organisms owing to the emphasis of long-term survey networks in the area, such as the extensive Multi-Agency Rocky Intertidal Network (MARINe for short). The vast majority of what we know about SSWD dynamics comes from studies along the Pacific Coast of North America using data produced through the combined effort of dedicated scientists, government agencies, and local citizens. The purpose of this investigation is to make use of this data monitoring system to investigate spatial and temoporal changes in sea star abundance along the west coast of the United States.
The Data
Below I explore the raw data I downloaded from a link I recieved after submitting a request for the data on the MARINe webpage.
library(readr)
raw_SS_count <- read.csv("~/Documents/GitHub/PaigeDuffin-Project/data/raw_data/seastarkat_size_count_totals_download.csv")
head(raw_SS_count, 3)## group_code site_code marine_site_name marine_sort_order latitude
## 1 UCLA ALEG Alegria 6420 34.46714
## 2 UCLA ALEG Alegria 6420 34.46714
## 3 UCLA ALEG Alegria 6420 34.46714
## longitude marine_common_season marine_season_code marine_common_year
## 1 -120.2774 85 SP02 2002
## 2 -120.2774 85 SP02 2002
## 3 -120.2774 85 SP02 2002
## season_sequence season_name target_assemblage method_code species_code
## 1 1 Spring sea_star IP PISOCH
## 2 1 Spring sea_star IP PISOCH
## 3 1 Spring sea_star IP PISOCH
## size_sort_order size_bin total mpa_designation mpa_region georegion
## 1 3 20 1 reference South Coast CA South
## 2 4 30 3 reference South Coast CA South
## 3 5 40 16 reference South Coast CA South
## bioregion state_province island last_updated
## 1 Government Point to Mexico California Mainland 2019-04-05 20:07:59
## 2 Government Point to Mexico California Mainland 2019-04-05 20:07:59
## 3 Government Point to Mexico California Mainland 2019-04-05 20:07:59levels(raw_SS_count$georegion)## [1] "AK" "CA Central"
## [3] "CA Channel Islands South" "CA North"
## [5] "CA North Central" "CA South"
## [7] "OR" "WA Olympic Coast"
## [9] "WA Salish Sea"levels(raw_SS_count$species_code)## [1] "EVATRO" "KATTUN" "PISOCH"Notably, it contains the following information (some parameters removed from this overview):
groupcode: The unique code for each monitoring groupmarine_site_name: The name of the site where the survey was conducted.latitudeandlongitude: (self explanatory)marine_season_code: A four-character code to identify the Sampling season. The first two characters indicate the season and the last 2 characters indicate the year.method_code: The method used for sampling the plot.size_bin: The size of the species being counted, binned to the nearest 5 or 10 millimeter.target_assemblage: The specific organism that this transect was created to monitor.total: Total number of individuals counted in a given size_binmpa_designation: Describes whether the referenced site is located within a Marine Protected Area (MPA) or is a reference site. If this field is blank, the referenced site is not located within an MPA, or is not a reference site for an existing MPA site.georegion/bioregion: geographic/biogeographic region in which site is locatedstate_province: The State or Province and Country where the referenced site is located.island: The name of the island where the referenced site is located. Sites not on islands are designated as mainland.
Additional Data Sources
I believe I have more data coming in to work with. I recieved an email today from a representative of the Multi-Agency Rocky Intertidal Network. This email said:
“You should have been able to download our sea star size/count data from the data download page that you were directed to after filling out the data request form. We also have associated disease category data, but because we only sample 1x per year at most sites, this additional info will not give you an idea of disease emergence or progression within a site or region. You might also be interested in our observation data set, which would be better for these types of questions. We can have this data to you within about a week.”
I’m also interested in environmental data such as sea surface temperature, dissolved oxygen, salinity, etc., as SSWD is thought to be exacerbated by abiotic stressors. However, I’m having a difficult time finding data to use that will pair reasonably well with the dates and locations of the raw_SS_count data I know I will be working with.
I’m very interested in this disease and environmental data, but don’t know what it will look like or if it will be feasible to use. Hopefully, I can incorporate this into my analysis, but for right now I am going to focus on outlining the questions I can answer with the data I have in hand.
Questions/Hypotheses to be addressed
State the research questions you plan to answer with this analysis
1. Are sea star populations changing over time? I anticipate that there will be a decline in sea star populations over time, but that the change will not occur uniformly across geographical regions and species. Species that are reported to have been affected by SSWD will experience the biggest declines.
2. Is species composition uniform across space? I expect to be able to define rough species ranges.
3. Is there evidence of sampling bias? Is there a relationship between important survey parameters (mainly, abundance data) and sampling method (ex. group conducting the survey, plot sampling method, the target assemblage). At the very least, I expect there will be bias based on the targeted organism of the survey; if a given survey aims to characterize the presence of one species of sea stars, it is likely that they will seek out that species and record higher counts.
Questions I may be able to answer with disease and/or environmental data:
1. Is there a relationship between sea star abundance over time and prevalence of SSWD? There absolutely should be a relationship, based on the body of evidence presented on studies that have focused on individual species.
2. Are there species of sea stars that are more affected by disease than others? The literature implies that this may be the case, but there also may be a bias based on species that researchers are more interested in.
3. Are there regions that are more affected by disease than others, regardless of species? The most recent epidemic appears to be so widespread that I suspect there won’t be any significant trends here, but it’s certainly worth checking.
4. Is there a relationship between sea star wasting disease and dissolved oxygen levels? A recent, but leading hypothesis in the field of SSWD suggests that it may be triggered specifically by low oxygen conditions. I’m very interested in exploring relationships between SSWD and DO levels- this is actually the main curiosity that lead me to focus in on this data set.
5. Is there a relationship between sea star wasting disease and [other water quality parameters I’m able to find data for]?
Part 2 Submission:
Recreating a Visualization using R Code
I have spent a very long time trying to figure out how to get all of the different components of the .Rproj to “talk” to each other (ie, write script in “processingscript.R”, have it save the processed script to the ./data/processed_data/ folder as an “RDS”) and I can’t figure it out. I looked at the example provided by Dr. Handel/his student in the second module, as well as some of the other students who have submitted already publically on the GitHub page.
In the interest of having something to turn in for this assignment on time, I am going to make this R Markdown file on my progress and turn it into a website, as this is something I do know how to do. I don’t want to complain, but I have to admit I am feeling very frusterated, overwhelmed, and defeated right now :(
Throughout this document, I’ve highlighted places where I’ve had an especially difficult time in font color=“red”>red font.
font color=“red”>Generally speaking, I do not know how to use the bibtex reference approach. I may be able to figure it out on my own, but truly am worried about how much of a time commitment this is all going to be. Are there any resources you can point me towards so I can learn how to use this system?
Processing Script
This is what I would include in my processingscript.R file.
processing script
this script loads the raw data, processes and cleans it and saves it as Rds file in the processed_data folder
load needed packages. make sure they are installed.
library(caret) # Classification and Regression Training## Warning: package 'caret' was built under R version 3.5.2## Warning: package 'ggplot2' was built under R version 3.5.2library(corrplot)
library(dplyr) # A Grammar of Data Manipulations## Warning: package 'dplyr' was built under R version 3.5.2library(forcats) # Tools for Working with Categorical Variables (Factors)## Warning: package 'forcats' was built under R version 3.5.2library(gdata) # Various R Programming Tools for Data Manipulation
library(ggplot2) # Create Elegant Data Visualizations Using the Grammar of Graphics
library(ggrepel) # Automatically Position Non-Overlapping Text Labels with ‘ggplot2’## Warning: package 'ggrepel' was built under R version 3.5.2library(ggridges) # Ridgeline Plots in ‘ggplot2’
library(ggthemes) # Extra Themes, Scales, and Geoms for ‘ggplot2’## Warning: package 'ggthemes' was built under R version 3.5.2library(gmodels) # Various R Programming Tools for Model Fitting
library(grid) # The Grid Graphics Package
library(gridExtra)
library(knitr) # A General-Purpose Package for Dynamic Report Generation in R## Warning: package 'knitr' was built under R version 3.5.2library(lmtest) # Testing Linear Regression Models## Warning: package 'lmtest' was built under R version 3.5.2library(plyr) # Tools for Splitting, Applying and Combining Data
library(purrr) # Functional Programming Tools## Warning: package 'purrr' was built under R version 3.5.2library(readr) # Read Rectangular Text Data
library(rpart)## Warning: package 'rpart' was built under R version 3.5.2library(rattle)
library(scales) # Scale Functions for Visualization
library(shiny) # Web Application Framework for R## Warning: package 'shiny' was built under R version 3.5.2library(skimr) # Compact and Flexible Summaries of Data## Warning: package 'skimr' was built under R version 3.5.2library(stringr) # Simple, Consistent Wrappers for Common String Operations
library(tibble) # Simple Data Frames## Warning: package 'tibble' was built under R version 3.5.2library(tidyr) # Easily Tidy Data with ‘spread()’ and ‘gather()’ Functions## Warning: package 'tidyr' was built under R version 3.5.2library(tidyverse) # Install and Load ‘Tidyverse’ Packages
library(usmap)## Warning: package 'usmap' was built under R version 3.5.2library(rworldmap)
library(ggmap)## Warning: package 'ggmap' was built under R version 3.5.2library(maps)
library(mapdata)
library(cowplot)## Warning: package 'cowplot' was built under R version 3.5.2library("googleway")
library("ggspatial")
library("rnaturalearth")
library("rnaturalearthdata")
theme_set(theme_bw())
library("sf")## Warning: package 'sf' was built under R version 3.5.2register_google(key = "AIzaSyB3n0SkJTXxrbUkF2bkBJs7XdleT_iJXsg", write = TRUE)
library(readxl)## Warning: package 'readxl' was built under R version 3.5.2library(dplyr)
library(readr)load data. path is relative to project directory.
# raw_SS_count <- read.csv("./data/raw_data/seastarkat_size_count_totals_download.csv") # already did that, above
head(raw_SS_count, 3)## group_code site_code marine_site_name marine_sort_order latitude
## 1 UCLA ALEG Alegria 6420 34.46714
## 2 UCLA ALEG Alegria 6420 34.46714
## 3 UCLA ALEG Alegria 6420 34.46714
## longitude marine_common_season marine_season_code marine_common_year
## 1 -120.2774 85 SP02 2002
## 2 -120.2774 85 SP02 2002
## 3 -120.2774 85 SP02 2002
## season_sequence season_name target_assemblage method_code species_code
## 1 1 Spring sea_star IP PISOCH
## 2 1 Spring sea_star IP PISOCH
## 3 1 Spring sea_star IP PISOCH
## size_sort_order size_bin total mpa_designation mpa_region georegion
## 1 3 20 1 reference South Coast CA South
## 2 4 30 3 reference South Coast CA South
## 3 5 40 16 reference South Coast CA South
## bioregion state_province island last_updated
## 1 Government Point to Mexico California Mainland 2019-04-05 20:07:59
## 2 Government Point to Mexico California Mainland 2019-04-05 20:07:59
## 3 Government Point to Mexico California Mainland 2019-04-05 20:07:59levels(raw_SS_count$georegion)## [1] "AK" "CA Central"
## [3] "CA Channel Islands South" "CA North"
## [5] "CA North Central" "CA South"
## [7] "OR" "WA Olympic Coast"
## [9] "WA Salish Sea"levels(raw_SS_count$species_code)## [1] "EVATRO" "KATTUN" "PISOCH"take a look at the data
dplyr::glimpse(raw_SS_count)## Observations: 13,165
## Variables: 24
## $ group_code <fct> UCLA, UCLA, UCLA, UCLA, UCLA, UCLA, UCLA, U…
## $ site_code <fct> ALEG, ALEG, ALEG, ALEG, ALEG, ALEG, ALEG, A…
## $ marine_site_name <fct> Alegria, Alegria, Alegria, Alegria, Alegria…
## $ marine_sort_order <int> 6420, 6420, 6420, 6420, 6420, 6420, 6420, 6…
## $ latitude <dbl> 34.46714, 34.46714, 34.46714, 34.46714, 34.…
## $ longitude <dbl> -120.2774, -120.2774, -120.2774, -120.2774,…
## $ marine_common_season <int> 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85,…
## $ marine_season_code <fct> SP02, SP02, SP02, SP02, SP02, SP02, SP02, S…
## $ marine_common_year <int> 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2…
## $ season_sequence <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ season_name <fct> Spring, Spring, Spring, Spring, Spring, Spr…
## $ target_assemblage <fct> sea_star, sea_star, sea_star, sea_star, sea…
## $ method_code <fct> IP, IP, IP, IP, IP, IP, IP, IP, IP, IP, IP,…
## $ species_code <fct> PISOCH, PISOCH, PISOCH, PISOCH, PISOCH, PIS…
## $ size_sort_order <int> 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 14, 15, 16…
## $ size_bin <int> 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 1…
## $ total <int> 1, 3, 16, 25, 27, 22, 11, 2, 2, 1, 1, 3, 6,…
## $ mpa_designation <fct> reference, reference, reference, reference,…
## $ mpa_region <fct> South Coast, South Coast, South Coast, Sout…
## $ georegion <fct> CA South, CA South, CA South, CA South, CA …
## $ bioregion <fct> Government Point to Mexico, Government Poin…
## $ state_province <fct> California, California, California, Califor…
## $ island <fct> Mainland, Mainland, Mainland, Mainland, Mai…
## $ last_updated <fct> 2019-04-05 20:07:59, 2019-04-05 20:07:59, 2…head(raw_SS_count, 3)## group_code site_code marine_site_name marine_sort_order latitude
## 1 UCLA ALEG Alegria 6420 34.46714
## 2 UCLA ALEG Alegria 6420 34.46714
## 3 UCLA ALEG Alegria 6420 34.46714
## longitude marine_common_season marine_season_code marine_common_year
## 1 -120.2774 85 SP02 2002
## 2 -120.2774 85 SP02 2002
## 3 -120.2774 85 SP02 2002
## season_sequence season_name target_assemblage method_code species_code
## 1 1 Spring sea_star IP PISOCH
## 2 1 Spring sea_star IP PISOCH
## 3 1 Spring sea_star IP PISOCH
## size_sort_order size_bin total mpa_designation mpa_region georegion
## 1 3 20 1 reference South Coast CA South
## 2 4 30 3 reference South Coast CA South
## 3 5 40 16 reference South Coast CA South
## bioregion state_province island last_updated
## 1 Government Point to Mexico California Mainland 2019-04-05 20:07:59
## 2 Government Point to Mexico California Mainland 2019-04-05 20:07:59
## 3 Government Point to Mexico California Mainland 2019-04-05 20:07:59levels(raw_SS_count$georegion)## [1] "AK" "CA Central"
## [3] "CA Channel Islands South" "CA North"
## [5] "CA North Central" "CA South"
## [7] "OR" "WA Olympic Coast"
## [9] "WA Salish Sea"levels(raw_SS_count$species_code)## [1] "EVATRO" "KATTUN" "PISOCH"change many of the variables from factors to characters, at least for now
SS_count <- raw_SS_count
SS_count$group_code <- as.character(SS_count$group_code)
SS_count$site_code <- as.character(SS_count$site_code)
SS_count$marine_site_name <- as.character(SS_count$marine_site_name)
SS_count$marine_season_code <- as.character(SS_count$marine_season_code)
SS_count$season_name <- as.character(SS_count$season_name)
SS_count$target_assemblage <- as.character(SS_count$target_assemblage)
SS_count$method_code <- as.character(SS_count$method_code)
SS_count$species_code <- as.character(SS_count$species_code)
SS_count$mpa_designation <- as.character(SS_count$mpa_designation)
SS_count$mpa_region <- as.character(SS_count$mpa_region)
SS_count$georegion <- as.character(SS_count$georegion)
SS_count$bioregion <- as.character(SS_count$bioregion)
SS_count$state_province <- as.character(SS_count$state_province)
SS_count$island <- as.character(SS_count$island)
SS_count$last_updated <- as.character(SS_count$last_updated)
glimpse(SS_count)## Observations: 13,165
## Variables: 24
## $ group_code <chr> "UCLA", "UCLA", "UCLA", "UCLA", "UCLA", "UC…
## $ site_code <chr> "ALEG", "ALEG", "ALEG", "ALEG", "ALEG", "AL…
## $ marine_site_name <chr> "Alegria", "Alegria", "Alegria", "Alegria",…
## $ marine_sort_order <int> 6420, 6420, 6420, 6420, 6420, 6420, 6420, 6…
## $ latitude <dbl> 34.46714, 34.46714, 34.46714, 34.46714, 34.…
## $ longitude <dbl> -120.2774, -120.2774, -120.2774, -120.2774,…
## $ marine_common_season <int> 85, 85, 85, 85, 85, 85, 85, 85, 85, 85, 85,…
## $ marine_season_code <chr> "SP02", "SP02", "SP02", "SP02", "SP02", "SP…
## $ marine_common_year <int> 2002, 2002, 2002, 2002, 2002, 2002, 2002, 2…
## $ season_sequence <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…
## $ season_name <chr> "Spring", "Spring", "Spring", "Spring", "Sp…
## $ target_assemblage <chr> "sea_star", "sea_star", "sea_star", "sea_st…
## $ method_code <chr> "IP", "IP", "IP", "IP", "IP", "IP", "IP", "…
## $ species_code <chr> "PISOCH", "PISOCH", "PISOCH", "PISOCH", "PI…
## $ size_sort_order <int> 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 14, 15, 16…
## $ size_bin <int> 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 1…
## $ total <int> 1, 3, 16, 25, 27, 22, 11, 2, 2, 1, 1, 3, 6,…
## $ mpa_designation <chr> "reference", "reference", "reference", "ref…
## $ mpa_region <chr> "South Coast", "South Coast", "South Coast"…
## $ georegion <chr> "CA South", "CA South", "CA South", "CA Sou…
## $ bioregion <chr> "Government Point to Mexico", "Government P…
## $ state_province <chr> "California", "California", "California", "…
## $ island <chr> "Mainland", "Mainland", "Mainland", "Mainla…
## $ last_updated <chr> "2019-04-05 20:07:59", "2019-04-05 20:07:59…skim(SS_count)## Skim summary statistics
## n obs: 13165
## n variables: 24
##
## ── Variable type:character ─────────────────────────────────────────────────────────
## variable missing complete n min max empty n_unique
## bioregion 0 13165 13165 13 33 0 6
## georegion 0 13165 13165 2 24 0 9
## group_code 0 13165 13165 2 6 0 9
## island 0 13165 13165 3 9 0 7
## last_updated 0 13165 13165 19 19 0 1
## marine_season_code 0 13165 13165 4 4 0 57
## marine_site_name 0 13165 13165 5 42 0 77
## method_code 0 13165 13165 2 4 0 4
## mpa_designation 0 13165 13165 2 9 0 5
## mpa_region 0 13165 13165 2 19 0 6
## season_name 0 13165 13165 4 6 0 4
## site_code 0 13165 13165 3 4 0 77
## species_code 0 13165 13165 6 6 0 3
## state_province 0 13165 13165 6 10 0 4
## target_assemblage 0 13165 13165 8 8 0 1
##
## ── Variable type:integer ───────────────────────────────────────────────────────────
## variable missing complete n mean sd p0 p25
## marine_common_season 0 13165 13165 116.82 19.23 78 101
## marine_common_year 0 13165 13165 2009.7 4.81 2000 2006
## marine_sort_order 0 13165 13165 5566.24 1125.15 1720 5020
## season_sequence 0 13165 13165 2.03 0.81 1 1
## size_bin 0 13165 13165 82.44 42.2 5 50
## size_sort_order 0 13165 13165 9.57 4.15 1 6
## total 0 13165 13165 9.92 16.43 1 2
## p50 p75 p100 hist
## 118 131 152 ▅▅▆▇▇▇▇▅
## 2010 2013 2018 ▃▅▅▇▇▆▆▆
## 6090 6370 7650 ▁▁▁▂▃▃▇▁
## 2 3 4 ▇▁▇▁▁▇▁▁
## 80 110 320 ▅▇▇▃▁▁▁▁
## 10 12 33 ▅▇▇▃▁▁▁▁
## 4 11 360 ▇▁▁▁▁▁▁▁
##
## ── Variable type:numeric ───────────────────────────────────────────────────────────
## variable missing complete n mean sd p0 p25 p50
## latitude 0 13165 13165 38.67 5.2 32.87 34.73 36.62
## longitude 0 13165 13165 -122.25 2.86 -135.38 -123.98 -121.94
## p75 p100 hist
## 41.65 57.05 ▇▆▃▂▁▁▁▁
## -120.62 -117.25 ▁▁▁▁▅▇▆▃save as RDS file
saveRDS(SS_count, file = “./data/processed_data/”)
font color=“red”>I don’t understand this step. After much trial and error, it continues to not work for me. Can you provide additional instruction on how to do this step? I keep recieving the error message: Error in gzfile(file, mode) : cannot open the connection. In addition: Warning message: In gzfile(file, mode) : cannot open compressed file ‘./data/processed_data/’, probable reason ‘Is a directory’, but when I write “./data/processed_data/processed_SS_count”, it gives me a different error
Exploring Study Sites
load necessary packages
library(caret) # Classification and Regression Training
library(corrplot)
library(dplyr) # A Grammar of Data Manipulations
library(forcats) # Tools for Working with Categorical Variables (Factors)
library(gdata) # Various R Programming Tools for Data Manipulation
library(ggplot2) # Create Elegant Data Visualizations Using the Grammar of Graphics
library(ggrepel) # Automatically Position Non-Overlapping Text Labels with ‘ggplot2’
library(ggridges) # Ridgeline Plots in ‘ggplot2’
library(ggthemes) # Extra Themes, Scales, and Geoms for ‘ggplot2’
library(gmodels) # Various R Programming Tools for Model Fitting
library(grid) # The Grid Graphics Package
library(gridExtra)
library(knitr) # A General-Purpose Package for Dynamic Report Generation in R
library(lmtest) # Testing Linear Regression Models
library(plyr) # Tools for Splitting, Applying and Combining Data
library(purrr) # Functional Programming Tools
library(readr) # Read Rectangular Text Data
library(rpart)
library(rattle)
library(scales) # Scale Functions for Visualization
library(shiny) # Web Application Framework for R
library(skimr) # Compact and Flexible Summaries of Data
library(stringr) # Simple, Consistent Wrappers for Common String Operations
library(tibble) # Simple Data Frames
library(tidyr) # Easily Tidy Data with ‘spread()’ and ‘gather()’ Functions
library(tidyverse) # Install and Load ‘Tidyverse’ Packages
library(usmap)
library(rworldmap)
library(ggmap)
library(maps)
library(mapdata)
library(cowplot)
library("googleway")
library("ggspatial")
library("rnaturalearth")
library("rnaturalearthdata")
theme_set(theme_bw())
library("sf")
register_google(key = "AIzaSyB3n0SkJTXxrbUkF2bkBJs7XdleT_iJXsg", write = TRUE)## Replacing old key (AIzaSyB3n0SkJTXxrbUkF2bkBJs7XdleT_iJXsg) with new key in /Users/Paige/.RenvironOverall histogram of counts per sample survey
SS_count %>% ggplot() +
geom_histogram(aes(total))## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Create new data frame that just lists each site once
unique_sites <- subset(raw_SS_count, !duplicated(raw_SS_count$site_code))
sites_coordinates <- unique_sites[c(2,5:6,15,22)]
sites_coordinates_contig <- filter(sites_coordinates, sites_coordinates$state_province != "Alaska")
CA_coordinates <- filter(sites_coordinates, sites_coordinates$state_province == "California")
OR_coordinates <- filter(sites_coordinates, sites_coordinates$state_province == "Oregon")
WA_coordinates <- filter(sites_coordinates, sites_coordinates$state_province == "Washington")
AK_coordinates <- filter(sites_coordinates, sites_coordinates$state_province == "Alaska")
CA_coordinates <- as.data.frame(CA_coordinates)
OR_coordinates <- as.data.frame(OR_coordinates)
WA_coordinates <- as.data.frame(WA_coordinates)
AK_coordinates <- as.data.frame(AK_coordinates)Plot sample sites
After messing around with several different mapping tools in R, I decided that trying to display the whole range (CA to AK) on one map would be really hard, and the Alaska sites are really just at one specific location, so I’m going to split Alaska off from the three states that are part of the contiguous US #### Plot AK sites
world <- ne_countries(scale = "medium", returnclass = "sf")
AK_zoom1 <- ggplot(data = world) +
geom_sf() +
geom_point(data= sites_coordinates,aes(x=sites_coordinates$longitude, y=sites_coordinates$latitude,), size = 2, alpha = 0.4) +
annotation_scale(location = "bl", width_hint = 0.3) +
coord_sf(xlim = c(-170, -128), ylim = c(54, 72), expand = FALSE) +
xlab("Longitude") +
ylab("Latitude") +
scale_x_continuous(breaks = c(-160, -150, -140))
AK_zoom2 <- ggplot(data = world) +
geom_sf() +
geom_point(data= sites_coordinates,aes(x=sites_coordinates$longitude, y=sites_coordinates$latitude,), size = 2, alpha = 0.4) +
annotation_scale(location = "bl", width_hint = 0.3) +
coord_sf(xlim = c(-144, -130), ylim = c(55, 62), expand = FALSE) +
xlab("Longitude") +
ylab("Latitude") +
scale_x_continuous(breaks = c(-142, -138, -134))
grid.arrange(AK_zoom1, AK_zoom2, ncol = 2)## Scale on map varies by more than 10%, scale bar may be inaccurate
## Scale on map varies by more than 10%, scale bar may be inaccurate
Plot WA, OR, and CA sites
State <- sites_coordinates_contig$state_province
ggplot(data = world) +
geom_sf() +
geom_jitter(data= sites_coordinates_contig,aes(x=sites_coordinates_contig$longitude, y=sites_coordinates_contig$latitude,col=State), size = 2, alpha = 0.6, width = 0.2) +
annotation_scale(location = "bl", width_hint = 0.3) +
coord_sf(xlim = c(-128, -115), ylim = c(32, 49.5), expand = FALSE) +
xlab("Longitude") +
ylab("Latitude")## Scale on map varies by more than 10%, scale bar may be inaccurate
Exploring Count Data by State and Year Surveyed
What does the structure of the data look like by state?
CA
CA <- SS_count %>% select(marine_site_name, marine_common_year, size_bin, total, state_province) %>%
filter(state_province == "California")
summary(CA)## marine_site_name marine_common_year size_bin total
## Length:10548 Min. :2000 Min. : 5.00 Min. : 1.000
## Class :character 1st Qu.:2005 1st Qu.: 50.00 1st Qu.: 2.000
## Mode :character Median :2009 Median : 80.00 Median : 5.000
## Mean :2009 Mean : 81.92 Mean : 9.607
## 3rd Qu.:2013 3rd Qu.:110.00 3rd Qu.: 11.000
## Max. :2018 Max. :230.00 Max. :360.000
## state_province
## Length:10548
## Class :character
## Mode :character
##
##
## OR
OR <- SS_count %>% select(marine_site_name, marine_common_year, size_bin, total, state_province) %>%
filter(state_province == "Oregon")
summary(OR)## marine_site_name marine_common_year size_bin total
## Length:1386 Min. :2000 Min. : 5.00 Min. : 1.00
## Class :character 1st Qu.:2005 1st Qu.: 50.00 1st Qu.: 2.00
## Mode :character Median :2010 Median : 80.00 Median : 7.00
## Mean :2010 Mean : 87.21 Mean : 15.28
## 3rd Qu.:2014 3rd Qu.:120.00 3rd Qu.: 18.00
## Max. :2018 Max. :230.00 Max. :212.00
## state_province
## Length:1386
## Class :character
## Mode :character
##
##
## WA
WA <- SS_count %>% select(marine_site_name, marine_common_year, size_bin, total, state_province) %>%
filter(state_province == "Washington")
summary(WA)## marine_site_name marine_common_year size_bin total
## Length:865 Min. :2008 Min. : 5.00 Min. : 1.000
## Class :character 1st Qu.:2011 1st Qu.: 50.00 1st Qu.: 1.000
## Mode :character Median :2014 Median : 90.00 Median : 2.000
## Mean :2014 Mean : 89.26 Mean : 6.776
## 3rd Qu.:2016 3rd Qu.:120.00 3rd Qu.: 5.000
## Max. :2018 Max. :320.00 Max. :133.000
## state_province
## Length:865
## Class :character
## Mode :character
##
##
## AK
AK <- SS_count %>% select(marine_site_name, marine_common_year, size_bin, total, state_province) %>%
filter(state_province == "Alaska")
summary(AK)## marine_site_name marine_common_year size_bin total
## Length:366 Min. :2011 Min. : 5.00 Min. : 1.00
## Class :character 1st Qu.:2013 1st Qu.: 40.00 1st Qu.: 1.00
## Mode :character Median :2015 Median : 60.00 Median : 3.00
## Mean :2015 Mean : 63.51 Mean : 6.03
## 3rd Qu.:2017 3rd Qu.: 87.50 3rd Qu.: 7.00
## Max. :2018 Max. :190.00 Max. :117.00
## state_province
## Length:366
## Class :character
## Mode :character
##
##
## [INSERT TABLE SUMMARIZING THIS INFORMATION]
What does sea star abundance over time look like?
SS_count %>% ggplot(aes(marine_common_year, total)) + geom_jitter(width = 0.15, alpha = 0.25)
Lets color the points by state
SS_count %>% ggplot(aes(marine_common_year, total, col=state_province)) + geom_jitter(width = 0.15, alpha = 0.25)
Woah, California is dominating this visualization.
Lets break that up into 4 parts:
SS_count %>%
ggplot(aes(marine_common_year, total, col=state_province)) + geom_jitter(width = 0.15, alpha = 0.25) +
facet_wrap(~state_province)
Washington and Alaska clearly weren’t surveyed for all years during which California and Oregon were.
What do these counts look like if we display them as log()?
SS_count %>%
ggplot(aes(marine_common_year, log(total), col=state_province)) + geom_jitter(width = 0.15, alpha = 0.25) +
facet_wrap(~state_province)
Abundance over time
SS_count %>% ggplot(aes(marine_common_year, log(total), group = marine_common_year)) + geom_boxplot()
So far I’ve been looking at this across all 3 species of sea stars surveyed. What does the last graph look like for each species, seperately?
Summary statistics for each species over time
PISOCH_only <- filter(SS_count, species_code == "PISOCH")
KATTUN_only <- filter(SS_count, species_code == "KATTUN")
EVATRO_only <- filter(SS_count, species_code == "EVATRO")
skim(PISOCH_only)## Skim summary statistics
## n obs: 11416
## n variables: 24
##
## ── Variable type:character ─────────────────────────────────────────────────────────
## variable missing complete n min max empty n_unique
## bioregion 0 11416 11416 13 33 0 6
## georegion 0 11416 11416 2 24 0 9
## group_code 0 11416 11416 2 6 0 9
## island 0 11416 11416 3 9 0 7
## last_updated 0 11416 11416 19 19 0 1
## marine_season_code 0 11416 11416 4 4 0 57
## marine_site_name 0 11416 11416 5 42 0 77
## method_code 0 11416 11416 2 4 0 4
## mpa_designation 0 11416 11416 2 9 0 5
## mpa_region 0 11416 11416 2 19 0 6
## season_name 0 11416 11416 4 6 0 4
## site_code 0 11416 11416 3 4 0 77
## species_code 0 11416 11416 6 6 0 1
## state_province 0 11416 11416 6 10 0 4
## target_assemblage 0 11416 11416 8 8 0 1
##
## ── Variable type:integer ───────────────────────────────────────────────────────────
## variable missing complete n mean sd p0 p25
## marine_common_season 0 11416 11416 114.4 19.14 78 99
## marine_common_year 0 11416 11416 2009.09 4.79 2000 2005
## marine_sort_order 0 11416 11416 5682.02 1048.93 1720 5050
## season_sequence 0 11416 11416 2.03 0.84 1 1
## size_bin 0 11416 11416 86.83 42.67 5 50
## size_sort_order 0 11416 11416 9.68 4.27 1 6
## total 0 11416 11416 10.85 17.41 1 2
## p50 p75 p100 hist
## 114 130 152 ▅▅▆▇▆▆▅▃
## 2009 2013 2018 ▃▅▅▇▆▆▃▅
## 6150 6390 7650 ▁▁▁▂▃▂▇▁
## 2 3 4 ▇▁▆▁▁▇▁▁
## 90 120 320 ▅▇▇▅▁▁▁▁
## 10 13 33 ▅▇▇▅▁▁▁▁
## 5 13 360 ▇▁▁▁▁▁▁▁
##
## ── Variable type:numeric ───────────────────────────────────────────────────────────
## variable missing complete n mean sd p0 p25 p50
## latitude 0 11416 11416 38.13 4.82 32.87 34.55 36.51
## longitude 0 11416 11416 -121.94 2.57 -135.38 -123.51 -121.91
## p75 p100 hist
## 40.34 57.05 ▇▅▂▂▁▁▁▁
## -120.61 -117.25 ▁▁▁▁▅▇▇▅skim(KATTUN_only)## Skim summary statistics
## n obs: 1676
## n variables: 24
##
## ── Variable type:character ─────────────────────────────────────────────────────────
## variable missing complete n min max empty n_unique
## bioregion 0 1676 1676 13 33 0 4
## georegion 0 1676 1676 2 16 0 7
## group_code 0 1676 1676 2 6 0 6
## island 0 1676 1676 3 9 0 5
## last_updated 0 1676 1676 19 19 0 1
## marine_season_code 0 1676 1676 4 4 0 27
## marine_site_name 0 1676 1676 5 42 0 46
## method_code 0 1676 1676 2 2 0 1
## mpa_designation 0 1676 1676 2 9 0 5
## mpa_region 0 1676 1676 2 19 0 5
## season_name 0 1676 1676 4 6 0 3
## site_code 0 1676 1676 3 4 0 46
## species_code 0 1676 1676 6 6 0 1
## state_province 0 1676 1676 6 10 0 4
## target_assemblage 0 1676 1676 8 8 0 1
##
## ── Variable type:integer ───────────────────────────────────────────────────────────
## variable missing complete n mean sd p0 p25 p50
## marine_common_season 0 1676 1676 132.29 10.05 114 123 133
## marine_common_year 0 1676 1676 2013.57 2.54 2009 2011 2014
## marine_sort_order 0 1676 1676 4911.32 1220.64 1720 4400 5050
## season_sequence 0 1676 1676 2.01 0.62 1 2 2
## size_bin 0 1676 1676 53.79 23.85 5 35 50
## size_sort_order 0 1676 1676 8.95 3.05 1 7 9
## total 0 1676 1676 3.9 3.26 1 1 3
## p75 p100 hist
## 141 150 ▆▇▅▆▇▇▅▆
## 2016 2018 ▇▇▅▆▇▇▆▇
## 6130 6310 ▂▁▁▁▂▇▁▆
## 2 3 ▂▁▁▇▁▁▁▂
## 70 130 ▃▃▇▅▇▂▁▁
## 11 17 ▂▂▃▇▇▅▁▁
## 5 19 ▇▂▂▁▁▁▁▁
##
## ── Variable type:numeric ───────────────────────────────────────────────────────────
## variable missing complete n mean sd p0 p25 p50
## latitude 0 1676 1676 41.73 5.77 35.45 36.56 40.34
## longitude 0 1676 1676 -124.07 3.33 -135.38 -124.14 -123.73
## p75 p100 hist
## 44.24 57.05 ▇▅▆▂▃▁▁▂
## -121.95 -120.95 ▁▁▁▁▁▁▇▆skim(EVATRO_only)## Skim summary statistics
## n obs: 73
## n variables: 24
##
## ── Variable type:character ─────────────────────────────────────────────────────────
## variable missing complete n min max empty n_unique
## bioregion 0 73 73 13 26 0 2
## georegion 0 73 73 2 13 0 2
## group_code 0 73 73 2 6 0 5
## island 0 73 73 3 9 0 6
## last_updated 0 73 73 19 19 0 1
## marine_season_code 0 73 73 4 4 0 10
## marine_site_name 0 73 73 9 42 0 10
## method_code 0 73 73 2 2 0 1
## mpa_designation 0 73 73 4 4 0 1
## mpa_region 0 73 73 4 4 0 1
## season_name 0 73 73 4 6 0 3
## site_code 0 73 73 3 4 0 10
## species_code 0 73 73 6 6 0 1
## state_province 0 73 73 6 10 0 2
## target_assemblage 0 73 73 8 8 0 1
##
## ── Variable type:integer ───────────────────────────────────────────────────────────
## variable missing complete n mean sd p0 p25 p50
## marine_common_season 0 73 73 140.75 9.19 114 134 142
## marine_common_year 0 73 73 2015.7 2.28 2009 2014 2016
## marine_sort_order 0 73 73 2496.03 833.39 1720 1720 1740
## season_sequence 0 73 73 1.96 0.26 1 2 2
## size_bin 0 73 73 54.11 28.41 5 30 50
## size_sort_order 0 73 73 6.4 2.87 1 4 6
## total 0 73 73 1.84 1.83 1 1 1
## p75 p100 hist
## 150 150 ▁▁▁▂▃▁▃▇
## 2018 2018 ▁▁▁▂▃▁▃▇
## 3290 3850 ▇▁▁▁▁▃▃▁
## 2 3 ▁▁▁▇▁▁▁▁
## 70 160 ▅▇▆▅▅▁▁▁
## 8 17 ▅▇▆▅▅▁▁▁
## 2 13 ▇▁▁▁▁▁▁▁
##
## ── Variable type:numeric ───────────────────────────────────────────────────────────
## variable missing complete n mean sd p0 p25 p50 p75
## latitude 0 73 73 52.86 4.5 47.58 48.52 56.99 57.05
## longitude 0 73 73 -129.42 6.4 -135.38 -135.35 -135.32 -122.55
## p100 hist
## 57.05 ▇▁▁▁▁▁▁▇
## -122.52 ▇▁▁▁▁▁▁▇Species abundance over time
PISOCH_only %>% ggplot(aes(marine_common_year, log(total), group = marine_common_year)) + geom_boxplot()
KATTUN_only %>% ggplot(aes(marine_common_year, log(total), group = marine_common_year)) + geom_boxplot()
EVATRO_only %>% ggplot(aes(marine_common_year, log(total), group = marine_common_year)) + geom_boxplot()
That last one looked like the log() function might not be applicable. Let’s see what it looks like without that:
EVATRO_only %>% ggplot(aes(marine_common_year, total, group = marine_common_year)) + geom_boxplot()
Characterizing sea star abundance in the last decade
The most recent and catastrophic sea star wasting disease outbreak hit the west coast of North America around 2013. How did sea star abundance look before the epidemic, how did it change during the epidemic, and have the populations rebounded?
last_decade <- filter(SS_count, marine_common_year >= "2010")
last_decade %>% ggplot(aes(marine_common_year, log(total), group = marine_common_year)) + geom_boxplot()
It looks like populations leveled back out at 2016 and remained relatively unchanged in 2017 and 2018.
Narrowing the time frame under analysis
It looks like theres a detectable decrease between 2013 and 2014, and that populations leveled back out at 2016 and remained relatively unchanged in 2017 and 2018. Lets narrow the time frame to only include sampling between 2013 and 2016.
last_decade$marine_common_year <- as.factor(last_decade$marine_common_year)
SSWD_timeframe <- filter(SS_count, marine_common_year >= "2013" & marine_common_year <= "2016")
SSWD_timeframe %>% ggplot(aes(marine_common_year, log(total), group = marine_common_year, col=marine_common_year)) + geom_boxplot()
By state
SSWD_timeframe %>% ggplot(aes(marine_common_year, log(total), group = marine_common_year, col=marine_common_year)) + geom_boxplot() + facet_wrap(~state_province, 1) + theme(legend.position = "none")
Fitting a linear model
SSWD_timeframe %>% ggplot(aes(marine_common_year, log(total), col=state_province)) + geom_jitter(width = 0.45, alpha = 0.25) +
facet_wrap(~state_province, 1) + geom_smooth(method = "lm", col="black") + theme(legend.position = "none")
Just California
SSWD_timeframe_CA <- filter(SSWD_timeframe, state_province == "California")
SSWD_timeframe_CA %>% ggplot(aes(marine_common_year, log(total), col=log(total))) + geom_jitter(width = 0.5, alpha = 0.25) + geom_smooth(method = "lm", col="red")
Doing the same as above, but separating by species, with a primary focus on Pisaster ochraceus (the organism I’m studing in my research)
SSWD_timeframe_PIS <- filter(PISOCH_only, marine_common_year >= "2013" & marine_common_year <= "2016")
SSWD_timeframe_PIS %>% ggplot(aes(marine_common_year, log(total), group = marine_common_year, col=marine_common_year)) + geom_boxplot() + facet_wrap(~state_province, 1) + theme(legend.position = "none")
SSWD_timeframe_KAT <- filter(KATTUN_only, marine_common_year >= "2013" & marine_common_year <= "2016")
SSWD_timeframe_KAT %>% ggplot(aes(marine_common_year, log(total), group = marine_common_year, col=marine_common_year)) + geom_boxplot() + facet_wrap(~state_province, 1) + theme(legend.position = "none")
SSWD_timeframe_EVA <- filter(EVATRO_only, marine_common_year >= "2013" & marine_common_year <= "2016")
SSWD_timeframe_EVA %>% ggplot(aes(marine_common_year, log(total), group = marine_common_year, col=marine_common_year)) + geom_boxplot() + facet_wrap(~state_province, 1) + theme(legend.position = "none")
OH! I didn’t catch that the reason why the EVA species had such few survey sightings was because its range appears to only be in Alaska. Speaking of which….
I should go back and plot each survey entry and color by species observed to get an idea of the ranges of each.
Species <- SSWD_timeframe$species_code
ggplot(data = world) +
geom_sf() +
geom_jitter(data= SSWD_timeframe,aes(x=SSWD_timeframe$longitude, y=SSWD_timeframe$latitude,col=Species), size = 2, alpha = 0.3, width = 0.2) +
annotation_scale(location = "bl", width_hint = 0.3) +
coord_sf(xlim = c(-128, -115), ylim = c(32, 49.5), expand = FALSE) +
xlab("Longitude") +
ylab("Latitude")## Scale on map varies by more than 10%, scale bar may be inaccurate
AK_zoom1_SSWD <- ggplot(data = world) +
geom_sf() +
geom_point(data= SSWD_timeframe,aes(x=SSWD_timeframe$longitude, y=SSWD_timeframe$latitude,col=Species,), size = 2, alpha = 0.4) +
annotation_scale(location = "bl", width_hint = 0.3) +
coord_sf(xlim = c(-140, -130), ylim = c(55, 62), expand = FALSE) +
xlab("Longitude") +
ylab("Latitude")
AK_zoom2_SSWD <- ggplot(data = world) +
geom_sf() +
geom_point(data= SSWD_timeframe,aes(x=SSWD_timeframe$longitude, y=SSWD_timeframe$latitude,col=Species,), size = 2, alpha = 0.4) +
annotation_scale(location = "bl", width_hint = 0.3) +
coord_sf(xlim = c(-136, -135), ylim = c(56.6, 57.4), expand = FALSE) +
xlab("Longitude") +
ylab("Latitude")
grid.arrange(AK_zoom1_SSWD, AK_zoom2_SSWD, ncol = 2)## Scale on map varies by more than 10%, scale bar may be inaccurate
Other components (unfinished)
Methods and Results
In most research papers, results and methods are separate. You can combine them here if you find it easier. You are also welcome to structure things such that those are separate sections.
Data aquisition
As applicable, explain where and how you got the data. If you directly import the data from an online source, you can combine this section with the next.
Data import and cleaning
Write code that reads in the file and cleans it so it’s ready for analysis. Since this will be fairly long code for most datasets, it might be a good idea to have it in one or several R scripts. If that is the case, explain here briefly what each file does. The files themselves should be commented well so everyone can follow along.
Univariate analysis
Use a combination of text/tables/figures to explore and describe your data. You should produce plots or tables or other summary quantities for most of your variables. You definitely need to do it for the important variables, i.e. if you have main exposure or outcome variables, those need to be explored. Depending on the total number of variables in your dataset, explore all or some of the others.
Bivariate analysis
Create plots or tables and compute simple statistics (e.g. t-tests, simple regression model with 1 predictor, etc.) to look for associations between your outcome(s) and each individual predictor variable
Full analysis
Use one or several suitable statistical/machine learning methods to analyze your data and to produce meaningful figures, tables, etc. This might again be code that is best placed in one or several separate R scripts that need to be well documented. You can then load the results produced by this code
Discussion
Summary and Interpretation
Summarize what you did, what you found and what it means.
Strengths and Limitations
Discuss what you perceive as strengths and limitations of your analysis.
Conclusions
What are the main take-home messages?
Include citations in your Rmd file using bibtex, the list of references will automatically be placed at the end